Genome-wide identification of A-to-I RNA editing events provides the functional implications in PDAC

Introduction RNA editing, a wide-acknowledged post-transcriptional mechanism, has been reported to be involved in the occurrence and development of cancer, especially the abnormal alteration of adenosine to inosine. However, fewer studies focus on pancreaticcancer. Therefore, we aimed to explore the possible linkages between altered RNA editing events and the development of PDAC. Method We characterized the global A-to-I RNA editing spectrum from RNA and matched whole-genome sequencing data of 41 primary PDAC and adjacent normal tissues. The following analyses were performed: different editing level and RNA expression analysis,pathway analysis, motif analysis, RNA secondary structure analysis, alternative splicing events analysis, and survival analysis.The RNA editing of single-cell RNA public sequencing data was also characterized. Result A large number of adaptive RNA editing events with significant differences in editing levels were identified, which are mainly regulated by ADAR1. Moreover, RNA editing in tumors has a higher editing level and more abundant editing sites in general. 140genes were screened out since they were identified with significantly different RNA editing events and were significantly different in expression level between tumor and matched normal samples. Further analysis showed a preference that in the tumor-specific group, they are mainly enriched in cancer-related signal pathways, while in the normal tissue-specific group, they are mainly enriched in pancreatic secretion. At the same time, we also found positively selected differentially edited sites in a series of cancer immune genes, including EGF, IGF1R, and PIK3CD. RNA editing might participate in pathogenisis of PDAC through regulating the alternative splicing and RNA secondary structure of important genesto further regulate gene expression and protein synthesis, including RAB27B and CERS4. Furthermore, single cell sequencing results showed that type2 ductal cells contributed the most to RNA editing events in tumors. Conclusion RNA editing is an epigenetic mechanism involved in the occurrence and development of pancreatic cancer, which has the potential to diagnose of PDAC and is closely related to the prognosis.


Introduction
RNA editing, one of the most common and abundant forms of posttranscriptional RNA modifications, has a vital impact on many biological processes (1), including missense codon changes (2), modulation of alteration splicing, RNA secondary structure, or modification of regulatory RNAs and their binding sites. It has been declared that RNA editing is related to many cancers, and abnormal changes in RNA editing level can induce cancers (3), including hepatocellular carcinoma (4), chronic lymphocytic leukemia (5), and ovarian cancer (6). Various phenomena and evidence have suggested that RNA editing has important physiological and pathological significance. For example, Masayuki Sakurai's team has found that abnormally regulated RNA editing events in the central nervous system played an important role in neurological development and brain function and were related to the pathogenesis of neurological and psychiatric disorders (7). Peng's team has found that A-to-I RNA editing events in cancer is a novel source of cancer protein heterogeneity, which can promote the protein diversity of cancer cells by altering the amino acid sequences (8). Therefore, RNA editing events in cancer should be paid more attention.
In humans, the most prevalent RNA editing event is the conversion of adenosine (A) on the double-stranded RNAs to inosine (I) through oxidative deamination reaction under the catalysis of adenosine deaminase ADARs (9). There are three types of ADAR proteins in mammals, ADAR1, ADAR2 (ADARB1), and ADAR3 (ADARB2), of which ADAR1 and ADAR2 are the unique mediators of A-to-I RNA editing and exist in most human tissues. Recent research has indicated that the abnormal expression of ADARs can trigger various diseases (2,10). ADAR1 and ADAR2 are reported to be involved in the cell proliferation activity and inflammatory response (11,12). In Okugawa's study, the abnormal activity of RNA editing in tumors is closely related to the overexpression of ADAR1 (13).
Various studies show that RNA editing participates in the pathogenesis of tumors and has the potential as a clinical indicator. Chigaev's study has systematically characterized the RNA editing genomic landscape of various cancers conducted through The Cancer Genome Atlas (TCGA) project or Genotype-Tissue Expression (GTEx) Project (14), which presents a difference in RNA editing level between tumor and normal samples. Higher levels of RNA editing were identified in most tumor samples (for example, head and neck cancer, breast cancer, and thyroid cancer). However, only a few coding RNA editing sites have been characterized (4), such as AZIN1, GABRA3, FLNB, SLC22A3, and COPA, which can affect tumor progression. Stably edited AZIN1 acts as an analog of ornithine decarboxylase (ODC) to prevent the degradation of ODC and cyclin D1, resulting in increased cell proliferation, metastasis potential and tumor initiation in colorectal cancer (15), esophageal squamous cell carcinoma (16), and hepatocellular carcinoma (17). In cervical cancer, once a tumor-suppressor BLCAP is edited, it loses its ability to interact with and inactivate STAT3, thereby increasing cell proliferation (18).
Pancreatic cancer is one of the malignant tumors with the highest fatality rate, gradually increased incidence, difficult early diagnosis, poor treatment effect, and strong tumor heterogeneity. However, the pathogenesis of pancreatic cancer is still unclear. Here, we aim to explore the pathogenic mechanism of pancreatic ductal adenocarcinoma (PDAC) from the perspective of RNA editing. We obtained a global profile on RNA editing events of PDAC based on the paired RNA and WGS data from 41 patients and found specific RNA editing events, genes, and signal pathways related to PDAC. A series of differential RNA editing events (DREs) was identified to be located at the important function genes, including pancreatic secretion, tumor immune-related genes, most of which also showed significant prognosis associations in PDAC. Thus, RNA editing might take a role in pathogenesis and provide implications for clinical outcome in PDAC.

Patients and samples
The patients diagnosed with PDAC were recruited from Changhai Hospital, Shanghai. All patients have signed a written informed consent form for collection and use of samples. A total of 41 pairs of tumors and their distal normal tissues were obtained from surgical specimens and frozen at -80°C immediately before nucleic acid extraction. All procedures have complied with the Code of Ethics of the World Medical Association (Declaration of Helsinki) and the guidelines of the institutional review committee of the Shanghai Institute of Nutrition and Health, Chinese Academy of Sciences.

DNA/RNA extraction and sequencing
Total RNA was isolated using the Qiagen RNeasy Kit (Qiagen, Germantown, MD) and treated with DNase I (New England Biolabs). Total DNA was extracted using the QIAamp DNA Mini Kit (Qiagen, Germantown, MD). The concentration was quantified by the NanoDrop One (Thermo Fisher Scientific, MA, USA), and the RNA quality was analyzed by the 2100 Agilent Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA). Libraries were generated using a standard library preparation kit (KAPA) according to the manufacturer's instructions. Libraries were quantified on the Qubit4 (Thermo Fisher Scientific, MA, USA) with the High Sensitivity DNA Assay (Thermo Fisher Scientific, MA, USA). Samples were equimolar pooled and sequenced on the Illumina NovaSeq (Illumina, San Diego, CA, USA) in 150-bp pairedend mode.

Data preprocessing
The quality of raw reads was assessed by FastQC (https://github. com/s-andrews/FastQC), and subsequently the low-quality data were filtered by fastp (https://github.com/OpenGene/fastp). The hg19/GRCh37 of the human genome was used as reference. RNAseq reads were aligned onto the reference genome using STAR ( h t t p s : / / g i t h u b . c o m / a l e x d o b i n / S T A R ) (parameters: -outSAMstrandField intronMotif -outSAMattributes Allr e a d F i l e s C o m m a n d z c a to u t F i l t e r T y p e B y S J o u to u t F i l t e r M u l t i m a p N m a x 1a l i g n S J o v e r h a n g M i n 8 -alignSJDBoverhangMin 1 -outFilterMismatchNmax 999 -outFilterMismatchNoverLmax 0.04 -alignIntronMin 20 -alignIntronMax 1000000 -alignMatesGapMax 1000000). WGS reads were aligned by BWA (https://github.com/lh3/bwa). All alignments were duplicated by Picard MarkDuplicates (http:// broadinstitute.github.io/picard/), sorted, and indexed using SAMtools (https://github.com/samtools/samtools). Computationally, RNA editing is identified as a single-nucleotide base change between DNA and RNA. The RNA and DNA variants in each sample were both detected by REDItools (19) (https://github.com/BioinfoUNIBA/ REDItools2) under unstranded strategy. Then, the sites containing DNA variants were removed and the rest were considered as RNA editing sites. Meanwhile, we also used SPRINT (20) (https://github. com/jumphone/SPRINT) to validate these RNA editing sites detected by REDItools.

Quantification and comparison of RNA editing levels
The number of A-to-I RNA editing sites is huge, but most of the sites exhibit low editing levels and low sample coverage. These editing events are often caused by non-specific overediting, whereas most of these editing sites are neutral or slightly harmful. This poses a huge challenge for detecting editing sites with high confidence. Therefore, we further obtain credible RNA editing sites by satisfying stringent requirements (base quality score >30, total reads for each site in each sample ≥10, 1 > editing level for each site in each sample ≥0.1, remove sites with DNA mutation and multiple RNA variants, remove SNP sites) ( Figure 1A).
RNA editing sites were classified into three groups, shared, normal-specific, and tumor-specific group, which implicated that RNA editing can not only produce new cancer-specific changes but also reedit existing editing sites in normal tissues to a certain extent. For the shared group, all editing sites were occurred in both tumor and normal samples and filtered by sample coverage over 10. The editing site was considered significantly different with a Fisher's exact p value less than 0.05 based on editing level. For two specific groups, the editing sites were occurred only in normal or tumor samples and filtered by sample coverage (≥4).

Gene ontology enrichment analysis
We used Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) molecular pathway enrichment analysis for the genes with RNA editing sites between tumor and normal tissues through R package "clusterProfiler" (http://www. bioconductor.org/packages/release/bioc/html/clusterProfiler.html). Statistical significance was defined as p.adjust < 0.05.

Motif search by MEME
The consensus motif in the MEME suite of motif-based sequence analysis tools (version 5.3.3, https://meme-suite.org/meme/) was identified using two motif widths of 30 and 50.

Construction of the protein-protein interaction network and module
We used the database of Search Tool for Retrieval of Interacting Genes (STRING) (ver. 11.0, http://www.string-db.org/), an online tool to analyze protein network characteristics and make protein network visualization, to evaluate the network of protein-protein interaction (PPI) and calculate K-means for clustering.

RNA secondary structure analysis
To predict the RNA secondary structure alteration derived from individual RNA editing events, we used the extended surrounding region (± 1000 nt around the editing site) to search for minimum free energy (MFE) and partition function using RNAfold (http://rna.tbi. univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). Meanwhile, MXfold2 (https://github.com/mxfold/mxfold2) was also applied to predict the RNA secondary structure of RNA sequence surrounding the editing sites (± 250 nt).
F ST calculation F ST analysis, as a genetic signature of positive selection, was performed using VCFtools.

Statistical analysis
Statistical analysis was performed by GraphPad Prism 9 (GraphPad Software, San Diego, CA, USA) and R packages. A value of p < 0.05 was considered statistically significant. Data were expressed as mean ± SEM.

Global properties of the inferred RNA editing sites in tumor and normal samples
To better understand the relationship between the RNA editing and PDAC, we respectively delineated the RNA editing profiles of normal and tumor samples and found that PDAC showed a severe A-to-I RNA editing imbalance. We analyzed the composition of high confident editing events and found that these events were enriched in A-to-I alternation (95.11%), which showed significant differences between tumor and normal samples (REDItools, paired t-test, p = 0.012) ( Figure 1B). The distribution of RNA editing events in chromosome was also analyzed, and a preference in chromosomes 17 and 19 was observed ( Figure 1C, Figure S1C).
We then assessed whether RNA editing sites showed differential editing levels between tumor and normal samples. A total of 1,657,409 and 654,485 unique A-to-I RNA editing sites were identified by REDItools and SPRINT, respectively, in which 504,659 editing sites were overlapped. The average number of A-to-I editing sites in tumor (average = 81504) was slightly higher than that in normal tissues (average = 70812) (paired t-test, p = 0.051) ( Figure 1D, Figure S1D). A percentage of 56.10% of the samples showed a hyper average editing level in tumors than those in normal tissues without difference overall (REDItools, paired t-test, p = 0.567) ( Figure 1E, Figure S1E). We also observed several genes that have been reported to change the editing level in cancer and even drive early tumor invasion and metastasis, including FLNB, SLC22A3, and AZIN1, which had significantly different editing levels in PDAC.
We further performed gene annotation and found that the number of editing sites per gene in tumor tissue was significantly higher than that in normal tissue (paired t-test, p = 0.002) ( Figure 1F). The distribution map showed that most genes have multiple editing sites ( Figure 1G). The types of RNA region showed an enrichment in intronic (80.99%), intergenic (9.29%), 3′UTR (3.07%), and exonic (4.06%) ( Figure 1H).

RNA editing showed different mechanisms in tumors and normal tissues
From a quantitative and qualitative perspective, the imbalance of RNA editing detected by REDItools was divided into three modes, shared group, tumor-specific group, and normal-specific group, based on the distribution characteristics of RNA editing events that occurred in tumor and normal samples. For A-to-I editing events, 13,042 and 7,297 events (belonging to 3,816 and 2,594 genes, respectively) occurred in the tumor-specific and normal-specific groups, respectively (Figure 2A), whereas 9,343 events (belonging to 2,812 genes) with significant differences occurred in the shared group ( Figure 2A). These editing sites among three groups were recognized as different RNA editing sites (DREs), in which 95.68% were identified by SPRINT. The level of DREs mainly ranged from 10% to 30% ( Figure 2B). In the shared group, there were more sites with high editing levels in tumor samples. Moreover, in the specific group, a low level of editing sites was predominant. Compared with the public REDIportal database (http://srv00.recas.ba.infn.it/atlas/), 58.91% of all identified events were known and 79.35%, 85.19%, and 96.36% were known in the tumor-specific group, normal-specific group, and shared group, respectively ( Figure 2C). The relative abundance of DREs showed a preference for chromosomes 17 and 19 ( Figures 2D, E). The distribution of DREs in each gene was consistent with the total result ( Figure 2F). Meanwhile, the types of RNA region of DREs showed a consistent result with total editing sites ( Figures S1G, H), and low enrichment of DREs was found in TSS (FigureS1I). Eight A-to-I DREs were identified as non-synonymous mutation, which might lead to the change in protein (Table S1). Furthermore, the DREs in the shared group displayed a hyper mean editing level and a higher editing event count in tumor than normal tissues (paired t-test, editing level: p < 0.001, editing event count: p = 0.023) (Figures 2H, I). In addition, the profile of total 29,682 A-to-I DREs could distinguish tumor and normal tissues well ( Figure 2G).
To investigate whether the genes with A-to-I DREs are involved in the physiological process of pancreas or tumorigenesis, we obtained 108 enriched pathways for total 5,319 genes with DREs by functional enrichment analysis. Pancreatic cancer, adherens junction, and MAPK signaling pathway were enriched in the tumor-specific group, whereas pancreatic secretion only existed in the normalspecific group (FigureS2A). In general, RNA editing events under different modes displayed different functional mechanisms.

ADAR1 regulates RNA editing events in PDAC
To determine whether ADAR family members have effects on RNA editing imbalance in PDAC, we compared the expression level of ADAR1 and its relation with the editing level and the amount of editing sites.
The expression of ADAR1 in tumor was significantly higher than that in normal tissues, rather than ADAR2 and ADAR3 (paired t-test, ADAR1: p < 0.0001, ADAR2: p = 0.357; ADAR3: p = 0.303) ( Figure 3A). Meanwhile, TCGA and GTEx data analysis showed that the expression profiles of the ADAR family were consistent with our results ( Figure  S3A). The A-to-I mean editing level was positively correlated with ADAR1 expression level, whereas it was lightly negative correlated with ADAR2 (spearman test, ADAR1: p < 0.0001, coef = 0.2174, ADAR2: p = 0.0285, coef = 0.059, ADAR3: p = 0.2689, coef = 0.01526) ( Figure 3B). The global association detected between ADAR1 expression and average editing level was also observed at each editing site ( Figure 3C). Meanwhile, only ADAR1 observed a significant association with prognosis (ADAR1, p = 0.042; ADAR2, p = 0.29; ADAR3, p = 0.066) ( Figure S3B). We then explored the triplet codon characteristics of A-to-I editing sites and found sequence preference, especially upstream G underrepresentation and downstream G excess, which was consistent with the known ADAR1 motif ( Figure 3D).
Then, we further classified samples with ADAR expression between paired tumor and normal tissues into ADAR1-upregulated and -downregulated groups to explore the function of ADAR1 on regulating RNA editing events. We found that the mean editing level was positively correlated with the amount of editing sites in the ADAR1 upregulated group. However, the opposite result was observed in the ADAR1-downregulated group ( Figure 3E).
Above all, these results suggest that the differential expression of ADAR1 manipulate the A-to-I imbalance in tumor samples in PDAC.

The genes with functional DREs
In order to explore the associations between RNA editing and gene expression in PDAC, differential gene expression analysis based on RNA-Seq was performed and 2,022 different expression genes (DEGs) were identified, including 712 upregulated and 1,310 downregulated genes in tumor ( Figure 4A). The DEGs were also found to distinguish tumor and normal tissues, but not better than the DREs did ( Figure 4B). Pathway analysis showed that downregulated DEGs were enriched in pancreatic secretion, PPAR signaling pathway, protein digestion and absorption, and cAMP signaling pathway, whereas upregulated DEGs were enriched in ECMreceptor interaction and IL-17 signaling pathway ( Figures 4C, D). Notably, in the pathogenesis of PDAC, the expression of tumorrelated genes changed greatly.
In order to detect whether the RNA editing events affect gene expression, we focus on the DEGs with DREs. A total of 140 genes were filtered out ( Figure 5A), which included 40 genes with 73 editing sites in the shared group, 78 genes with 370 A-to-I editing sites in the tumor-specific group, and 61 genes with 266 A-to-I editing sites in the normal-specific group. These differentially edited sites accompanied by different expressions were named eDREs, and these genes were named as edDEGs. The genes in the shared group were enriched in pancreatic secretion; focal adhesion; protein digestion and absorption; regulation of actin cytoskeleton; alanine, aspartate, and glutamate metabolism; and ECM-receptor interaction ( Figure 5B). For genes in the tumor-specific group, they were enriched in ECM-receptor interaction, focal adhesion, central carbon metabolism in cancer, regulation of actin cytoskeleton, and others. For genes in the normal-specific group, KEGG pathways were enriched in the PPAR signaling pathway; alanine, aspartate, and glutamate metabolism; pancreatic secretion; and protein digestion and absorption. Furthermore, only the gene set in the tumor-specific group observed a significant association with prognosis (total 140 edDEGs: p = 0.1, shared group: p = 0.022, normal-specific group: p = 0.13, tumorspecific group, p = 0.001) ( Figure 5C). These results suggested that RNA editing might take a role in the pathogenesis of PDAC.
Then, the expression levels of the 140 edDEGs were validated using TCGA-PAAD datasets (http://gepia2.cancer-pku.cn/#index). A set of 93 edDEGs were filtered out to showed similar expression levels with our results and showed a significant association with prognosis (p = 0.013) ( Figure 5D). We then analyzed the association between the level of RNA editing and the level of gene expression and found that it could be either positive or negative. For example, the level of the RNA editing event in chromosome 4: 110846540 on EGF showed a significantly negative correlation with its expression level (p = 0.0122) ( Figure 5E).
Among them, 41 edDEGs were found to be significantly associated with prognosis in TCGA-PAAD dataset, the combination of which also showed a significant association with prognosis (p < 0.001) ( Figure 5F). Moreover, the PPI network of 41 edDEGs was divided into three clusters, in which ITGA2, ITGA3, ITGB6, and LAMC2, MUC4, and B3GNT3 were the nodes with most connectivity ( Figure 5G). These genes were reported to be related with the pathogenesis of PDAC. ITGA2 can promote cell migration, invasion, metastasis, and chemoresistance in PDAC through ECM remodeling (22). LAMC2 can promote cancer progression and gemcitabine resistance through modulation of EMT and ATPbinding cassette transporters in PDAC (23). ITGA3 has been reported with an upregulated expression in PDAC (24). MUC4promoted neural invasion is mediated by the axon guidance factor Netrin-1 in PDAC (25). In summary, RNA editing indeed affects the expression of important functional genes, suggesting an important role of RNA editing in the pathogenesis of PDAC.

Functional eDREs affect RNA stability
It was reported that RNA editing could affect RNA stability by altering the RNA secondary structure and alternative splicing (26). Thus, we analyzed the RNA sequences and the corresponding secondary structure close to the editing site. We predicted the potential motif surrounding 177 eDREs derived from the 41 edDEGs by MEME and found the preference of GA-rich and CArich regions in sequences ( Figure 5H). Furthermore, we also found that 42.67% of these eDREs might impact the RNA secondary structure based on RNAfold analysis based on minimum free energy (MFE) and partition function (Table S2), which in turn might affect RNA stability and function. For example, the predicted RNA structure was changed from stem to loop at a GA-rich region due to the nearby RNA editing event in chromosome 11: 117986470 on TMPRSS4 (Figure 5I), which has been reported as an independent risk factor in PDAC to promote cell proliferation and inhibit apoptosis through activating ERK1/2 Signaling pathway (27).
A total of 3,647 significantly different alternative splicing events were identified between tumor and normal, of which 41.98% contained at least one RNA editing site ( Figure 5J, Table S3). These differentially edited sites accompanied by different alternative splicing events were named sDREs, and these genes were named edDSGs. Major sDREs participated in skipped exons. Among them, LAMC2, RAB27B, PHYHD1, MYH11, and CERS4 were also occurred in 140 edDEGs with co-changes in expression and editing. The alternative splicing events of LAMC2, RAB27B, MYH11, and CERS4 showed that once an RNA editing event occurred on the upstream intron, at least one downstream exon was shipped ( Figure 5L). Exon 3 of CERS4 is a protein-coding region, and its loss because of A-to-I editing in intron2 would result in truncated protein, thereby affecting its protein function. Two kinds of sDREs in intron 2 of RAB27B lead to an alteration of the RNA secondary structure and also generated an abnormal transcript ( Figure 5L). PHYHD1 displayed alternative 5′ splicing sites due to A-to-I editing in intron 10 leading to a truncated exon 10 ( Figure 5L). We further observed the editing levels among various alternative splicing types and found that hyper editing events accounted for a large proportion of each component in the tumor ( Figure 5K). These suggested that aberrant RNA editing in tumors may lead to aberrant variable splicing through an unstable RNA secondary structure. These results suggest that RNA editing regulates gene expression and protein synthesis by affecting alternative splicing and the RNA secondary structure.

Positive selection pressure on the RNA editing events
It was reported that positive selection is one of the mechanisms driving RNA editing (28,29). In order to determine whether the RNA editing events suffer from positive selection in the tumorigenesis, we performed calculation of F ST for each DREs. A percentage of 5.48% of DREs belonging to 1,005 genes showed the signal of positive selection ( Figure 6A). These genes were enriched in regulation of actin cytoskeleton, focal adhesion, and adherens junction ( Figure 6B), and 48 genes were immune-related genes. Furthermore, we screened out the positive selection signals of all immune genes (Table S4). Chromosomes 19 and 17 had relatively more selection signals ( Figure 6C), which was consistent with the enrichment of editing events in chromosomes 19 and 17. Meanwhile, an interesting result was found in which the positively selected DREs were mainly located in the normal-and tumor-specific groups. The hot genes with the most positive signal events in the normal-specific group were RBPJ, PFKFB3, CAMK1D, ACACB, and WWOX, whereas the hot genes in the tumor-specific group were SLC35F3, PCDH7, NF1, and RAB27B. Among them, RBPJ was an important transcriptional regulator in the Notch signaling pathway (30). PFKFB3 was required for cell-cycle progression and prevention of apoptosis (31).
RAB27B expression is considered an independent prognostic marker for PDAC (32).
In summary, positive selection indeed has an effect on the RNA editing in PDAC and has a preference for the RNA editing sites in tumor immune-related genes.

RNA editing events has a cell type bias
We downloaded a single-cell RNA-seq data set (PRJCA001063) of pancreatic cancer containing 24 tumors and 11 normal pancreases and identified 10 major clusters including type 1 ductal, type 2 ductal, acinar, endocrine, endothelial, fibroblast, stellate, macrophage, T, and B cells, through principal component analysis based on gene expression. We then extracted the RNA editing events for each cluster. In general, the significant difference in editing level occurred in type 2 ductal, macrophage, and T cells, whereas the significant difference in editing sites occurred in type 1 ductal, type 2 ductal, macrophage, stellate, and T cells ( Figure 7A). Among them, type 2 ductal cells greatly contributed to RNA editing events in pancreatic cancer tissues, whereas type 1 ductal cells had the largest number of RNA editing events in normal tissues. We then identified the DREs and their located genes and found that type 2 ductal, endothelial, fibroblast, stellate, macrophage, and T cells were dominated in tumor-specific editing sites, whereas normalspecific editing sites were dominated type 1 ductal cells ( Figure 7B). The shared group showed a significantly higher editing level in type1 ductal cells and a significantly lower editing level in fibroblast, stellate, and macrophage in tumors, whereas in the specific group, a significantly higher editing level occurred in acinar, B, type 2 ductal, endocrine, macrophage, and T cells and a significantly lower editing level was displayed in stellate cells ( Figure 7C). The Pearson analysis showed that it could be distinguished into normal and tumor based on the level of DREs ( Figure 7D). Then, we compared the distribution of genes with RNA editing sites between RNA-seq of 41 paired PDAC samples and single-cell sequencing data of 35 PDAC ( Figure 7E). There were 187 overlapped observed only in the tumor-specific group between RNAseq and single-cell sequencing, of which 72.06% were detected in type 2 ductal cells.
In conclusion, the level and amount of RNA editing events were different among various cell types of PDAC, and type 2 ductal cells contributed the most to RNA editing in tumors.

Discussion
A-to-I RNA editing is widespread in the human transcriptome, and the genetic variation generated by it can expand the diversity and complexity of the transcriptome. Recent studies have proved that RNA editing participates in the pathogenesis of cancer. In this study, we first focus on the pathogenic mechanism of PDAC from the perspective of RNA editing. Using matched genomic and transcriptomic data from 41 patients, we have delineated the comprehensive A-to-I RNA editing landscape and studied the distribution features of DREs and the associations with gene expression. We have observed that the spectrum of DREs can distinguish tumor and normal tissues well. The analysis of DREs and related genes suggests that RNA editing events are mainly regulated by ADAR1 and are involved in PDAC tumorigenesis. We have also found some RNA editing events that are functional and positively selected in tumor, and the genes with these editing events are closely related to the prognosis.
We have verified that RNA editing is mainly regulated by ADAR1 in PDAC. First, the vast majority of editing sites found fall in the ALU repeat sequence (Figure S1F), which is consistent with the reported preference for editing regions of ADAR1 (33). Second, ADAR1 is significantly highly expressed in tumors and is significantly positively correlated with the RNA editing level in PDAC, which is also consistent with the reported results. For example, Leng Han has revealed that the diversity of RNA editing events in tumor samples has the best correlation with the overall expression level of ADAR1 (1). A recent study has determined the frequency of RNA editing events in 17 cancer types in The Cancer Genome Atlas (TCGA) database 3 (34), in which the level of global RNA editing is positively correlated with ADAR1. Third, the triplet codon sequences centered on editing site A have a preference with 5′ U and 3′ G, which corresponds to the known motif of ADAR1 (35). Furthermore, several recent studies have emphasized the role of ADAR1 in cancer development. Chen found that ADAR1 edits multiple sites in the YXXQ motif of BLCAP (a tumor suppressor for bladder cancer), causing it to lose its inhibition of STAT3 activation, thereby promoting tumorigenesis in cervical cancer (CC) (18). Salameh found that ADAR1-mediated editing of prostate cancer antigen 3 (PCA3) increased its stability and expression in prostate cancer (PC), further promoting tumorigenesis (36). Meanwhile, we have also observed that some well-acknowledged pancreatic cancer genes contains DREs, such as BRCA2 (37), ATM (38), and SMAD4 (39). In summary, RNA editing is mainly regulated by ADAR1 to be involved in the occurrence and development of PDAC.
RNA editing events have been reported to suffer under selection pressure. Zhang Rui's team has found that RNA-edited SNP sites are highly enriched near autoimmune disease-related sites represented by Crohn's disease, asthma, and allergic dermatitis and are subject to balanced selection (40). Duan's group has reported that there are a large number of non-synonymous RNA editing sites (Nonsyn) that change amino acids in Drosophila, which show adaptive signals and are subject to positive natural selection (28). We boldly speculate that natural selection has an effect on the process of ADAR1-regulated RNA editing events in PDAC. Interestingly, some positively selected RNA editing sites are located in immune genes. For example, in this study, two tumor-specific editing sites in PPARG, chr3:12404457 and chr3:12432364, were identified with positive selection. PPARG has been reported to be expressed in various tumor cells and shows significant association with prognosis in PDAC (p = 0.016). The positively selected RNA editing sites in cancer tissues are likely to be produced by self-regulation of constantly adapting to the microenvironment and changing the microenvironment during the survival and reproduction of cancer cells. In a word, the results further indicate that RNA editing may be driven by nature selection to play a role in pathogenesis of PDAC.
RNA editing may be involved in pathogenesis of PDAC by posttranscriptional regulation such as variable splicing and RNA secondary structure change. Sze Jing Tang et al. proposed a novel mechanism in which ADAR1-mediated editing can affect alternative splicing, which in turn affects protein binding (41). They found that ADAR1 specifically edits GA-rich ISS at intron 8 of CCDC15, leading to recruitment of SRSF7 to the edited region and repression of exon 9 inclusion. Our study also found some events, which might affect posttranscriptional regulation. A series of RNA editing sites that affect the level of RNA transcription and alter RNA secondary structure (Table  S2) have been identified. For example, the RNA editing at chromosome 11: 117986470 in TMPRSS4 can change the RNA secondary structure, and the expression level of TMPRSS4 is significantly correlated with the editing level. These results show that RNA editing events can take a role in pathogenesis of PDAC by changing the stability of the secondary structure and the expression levels of key genes. However, experiments are needed for further validation.
In addition, our results suggest that RNA editing might act as a novel indicator for the diagnosis and prognosis of pancreatic cancer to some extent. We have observed that the RNA editing profile can more clearly distinguish tumors and normal samples, implying that RNA editing has potential for diagnosis of PDAC. The DEGs related to prognosis of PDAC (Table S2) suggest that RNA editing has the potential to judge the prognosis of PDAC. Among these genes, 16 genes have been reported to be associated with PDAC, including TMPRSS4 (42), ITGA2, ITGA3 (43), GPRC5A (1), LAMC2 (23), ARNTL2 (44), RAB27B (32), and PADI1 (45). RAB27B, a member of the Rab family GTPases involved in vesicle trafficking, has been found to be involved in PDAC invasion in a few studies. The exosome secretion pathway regulated by RAB27B is considered to be a novel therapeutic target for PDAC (32). In our study, we found that RAB27B can change the secondary structure of mRNA through RNA editing and then generate the transcript missing exon 2, which eventually leads to the structural change of protein. The three genes, VSIG1, IGFL2, and PLS1, have not been reported to be related with PDAC. VSIG1, a cell adhesion protein of the immunoglobulin superfamily, is preferentially expressed in gastric, testicular, esophageal, and ovarian cancers (46). VSIG1 has been reported to be related to the metastatic behavior of various colon cancer cell types (47), and nuclear positivity of VSIG1 has been observed in all cases of distant metastasis of gastrointestinal stromal tumors (48). PLS1 has been reported to promote metastasis of colorectal cancer through the IQGAP1/Rac1/ERK pathway (49). IGFL2 has been identified as a member of 12 marker panel of cancer-associated fibroblasts associated with the progression of hepatocellular carcinoma (50). The roles of these four genes with significantly differential expression and RNA editing levels in pancreatic cancer deserve further study.

Conclusions
These findings further expand our understanding of the role of RNA editing events in the occurrence and development of tumors. The specific occurrence of events in tumors may also have the potential to serve as important markers for PDAC diagnosis and prognosis. However, the understanding of the function of RNA editing is not comprehensive enough, and what kind of RNA editing events in which PDAC cells play a critical role in the occurrence and development of PDAC still needs a large amount of data for further research.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, GSE172356, https://db.cngb.org/, CNP0003871.

Ethics statement
The studies involving human participants were reviewed and approved by Naval Medical University. The patients/participants provided their written informed consent to participate in this study.